 
Module RobMean

  ! robust estimator of arithmetical mean by
  ! Hogg in Launer, Wilkinson: Robustness in Statistics

  ! print debug informations?
  logical, parameter, private :: print = .false.

contains

!-------------------------------------------------------------------

  subroutine MADest(n,x,median,mad)

! estimate median and median of absoluted deviations

    use QMedian

    integer :: n
    real :: x(:),median,mad,xx(n)

    xx = x(1:n)
    call qmed(n,xx,n/2,median)
    xx = abs(x(1:n) - median)
    call qmed(n,xx,n/2+1,mad)

    if( print ) &
    write(*,*) "median, median of absoluted deviations : ",median, mad

  end subroutine MADest

!-----------------------------------------------------------------------

  subroutine am1(n,x,t,dt,psi,dpsi,maxit,istat)

    ! Newton's iterations

    interface
       function psi(x)
         real :: psi,x
       end function psi
       function dpsi(x)
         real :: dpsi,x
       end function dpsi
    end interface

    integer :: n,i,it,maxit,istat
    real :: x(:)
    real :: t,dt
    real :: d,r,s,s2,sum1,sum2,sum3,rp

    istat = 0; dt = -1.0

    if( n < 1 ) then                     ! a few data
       istat = 1
       t = 0.0
    endif

    if( n == 1 ) then               ! only one point, but correct case
       t = x(1)
       dt = 0.0
       return
    endif

    ! initial values    
    call MADest(n,x,t,s)
    s = s/0.6745
    s2 = s**2

    if( abs(s) < epsilon(s) )then    ! identical points on input
       t = sum(x(1:n))/n
!       write(*,*) x(1:n),t,s
       dt = 0.0
       return
    endif

    do it = 1,maxit

       ! corrector's estimation 
       sum1 = 0.0; sum2 = 0.0; sum3= 0.0;
       do i = 1,n
          r = (x(i) - t)/s
          rp = psi(r)
          sum1 = sum1 + rp
          sum2 = sum2 + dpsi(r)
          sum3 = sum3 + rp**2
       enddo
       if( abs(sum2) < epsilon(sum2) )then
          istat = 2               ! ?
          return
       endif
!       write(*,*) sum1,sum2,sum3
       d = s*sum1/sum2
       t = t + d
       dt = s2*n/(n-1)*sum3/sum2**2

       ! exit of iterations: 
       if( it > 2 .and. (d**2 < 0.01*dt .or. abs(d) < 10.0*epsilon(d)) ) exit

       if( print ) &
            write(*,*) "prumer, prirustek, chyba: ",t,d,sqrt(dt),n
    enddo

    ! deviation's estimation
    sum1 = 0.0; sum3 = 0.0
    do i = 1, n
       r = (x(i) - t)/s
       sum1 = sum1 +  psi(r)**2
       sum3 = sum3 + dpsi(r)
    enddo
    if( sum3 == 0.0 )then
       istat = 3                ! ?
       return
    endif
    dt = s2*n/(n-1)*sum1/sum3**2
    dt = sqrt(dt)

  end subroutine am1

!------------------------------------------------------------------------

  subroutine am2(n,x,t,dt,psi,maxit,istat)

    ! H-method

    interface
       function psi(x)
         real :: psi,x
       end function psi
    end interface

    integer :: n,i,it,maxit,istat
    real :: x(:)
    real :: t,dt
    real :: d,s,sum1,sum2,r,s2

    istat = 0; dt = -1.0;

    ! initial values    
    call MADest(n,x,t,s)
    s = s/0.6745
    s2 = s**2

    if( s == 0.0 .or. n <= 1) then
       istat = 1                    ! a few values or all the same
       return
    endif

    do it = 1,maxit
       sum1 = 0.0
       do i = 1,n
          sum1 = sum1 + psi((x(i) - t)/s)
       enddo
       d = s*sum1/n
       t = t + d

       ! exit of iterations: 
       if( abs(d) < 0.1*dt .or. abs(d) < 10.0*epsilon(d) ) exit

       ! deviation estimator
       sum1 = 0.0; sum2 = 0.0
       do i = 1, n
          r = abs(psi((x(i) - t)/s))
          sum1 = sum1 + r*(x(i) - t)**2
          sum2 = sum2 + r
       enddo
       dt = sqrt( sum1/sum2/(n-1) )   

       if( print ) &
       write(*,*) "prumer, pocet hodnot, prirustek: ",t,d,dt
    enddo
 
  end subroutine am2

!--------------------------------------------------------------------

  subroutine am3(n,x,t,dt,psi,maxit)

    ! least-square with flow weigths
    
    interface
       function psi(x)
         real :: psi,x
       end function psi
    end interface

    integer :: n,i,it,maxit
    real :: x(:)
    real :: t,dt
    real :: d,s,sum1,sum2,r,psi1

    ! initial values    
    call MADest(n,x,t,s)    
    s = s/0.6745

    do it = 1,maxit
       sum1 = 0.0; sum2 = 0.0
       do i = 1,n
          r = (x(i) - t)/s
          psi1 = psi(r)
          sum1 = sum1 + psi1
          if( r /= 0.0 )then
             sum2 = sum2 + psi1/r
          endif
       enddo
       d = s*sum1/sum2
       t = t + d

       ! exit of iterations: 
       if( abs(d) < 0.1*dt .or. abs(d) < 10.0*epsilon(d) ) exit

       ! deviation estimator
       sum1 = 0.0; sum2 = 0.0
       do i = 1, n
          r = abs(psi((x(i) - t)/s))
          sum1 = sum1 + r*(x(i) - t)**2
          sum2 = sum2 + r
       enddo
       dt = sqrt( sum1/sum2/(n-1) )

       if( print ) &
       write(*,*) "prumer, pocet hodnot, prirustek: ",t,d,dt
    enddo
 
end subroutine

end module



